getwd()
setwd(getwd())
#setwd("D:\thasegawa\ExtremeEvent\R")
library(reshape2)
library(ggplot2)
library(plyr)
library(grid)
library(scales)
library(gdxrrw)
library(magrittr)

colnames <- c("VAR","Year","Region","SSP","GCM","RCP","CO2fert","SYR","N","Crop","Value")
colnamesNocc <- c("VAR","Year","Region","SSP","Crop","Value")
linepalette <- c("#4DAF4A","#377EB8","#E41A1C","#FF7F00","#984EA3","#FFFF33","#A65628","#F781BF","#8DD3C7","#FB8072","#80B1D3","#FDB462","#B3DE69","#FCCDE5","#D9D9D9","#BC80BD","#CCEBC5","#FFED6F")

regionlabel <- function(y){
  levels(y$Region)[levels(y$Region)=="WLD"] <- "World"
  levels(y$Region)[levels(y$Region)=="USA"] <- "United States"
  levels(y$Region)[levels(y$Region)=="XE25"] <- "EU25"
  levels(y$Region)[levels(y$Region)=="XER"] <- "Rest of Europe"
  levels(y$Region)[levels(y$Region)=="TUR"] <- "Turkey"
  levels(y$Region)[levels(y$Region)=="XOC"] <- "Oceania"
  levels(y$Region)[levels(y$Region)=="CHN"] <- "China"
  levels(y$Region)[levels(y$Region)=="IND"] <- "India"
  levels(y$Region)[levels(y$Region)=="JPN"] <- "Japan"
  levels(y$Region)[levels(y$Region)=="XSE"] <- "Southeast Asia"
  levels(y$Region)[levels(y$Region)=="XSA"] <- "Rest of South Asia"
  levels(y$Region)[levels(y$Region)=="CAN"] <- "Canada"
  levels(y$Region)[levels(y$Region)=="BRA"] <- "Brazil"
  levels(y$Region)[levels(y$Region)=="XLM"] <- "Rest of South America"
  levels(y$Region)[levels(y$Region)=="CIS"] <- "CIS"
  levels(y$Region)[levels(y$Region)=="XME"] <- "Middle East"
  levels(y$Region)[levels(y$Region)=="XNF"] <- "North Africa"
  levels(y$Region)[levels(y$Region)=="XAF"] <- "Sub-Sa Africa"
z <- y
return(z)
}

regionlabelagg <- function(y){
  levels(y$Region)[levels(y$Region)=="WLD"] <- "World"
  levels(y$Region)[levels(y$Region)=="US"] <- "United States"
  levels(y$Region)[levels(y$Region)=="LAM"] <- "Latin America"
  levels(y$Region)[levels(y$Region)=="EU"] <- "EU"
  levels(y$Region)[levels(y$Region)=="CHN"] <- "China"
  levels(y$Region)[levels(y$Region)=="SEAsia"] <- "Southeast Asia"
  levels(y$Region)[levels(y$Region)=="SAsia"] <- "South Asia"
  levels(y$Region)[levels(y$Region)=="Afr_ME"] <- "Africa & Middle East"
  levels(y$Region)[levels(y$Region)=="FSU"] <- "CIS"
  levels(y$Region)[levels(y$Region)=="Other"] <- "Other"
  z <- y
  return(z)
}
croplabel <- function(y){
  levels(y$Crop)[levels(y$Crop)=="PDR"] <- "rice"
  levels(y$Crop)[levels(y$Crop)=="WHT"] <- "wheat"
  levels(y$Crop)[levels(y$Crop)=="GRO"] <- "other grains"
  levels(y$Crop)[levels(y$Crop)=="OSD"] <- "oil crops"
  z <- y
  return(z)
}
fsize <- 13
MyThemeLine <- theme_grey() +
  theme(
    panel.border=element_rect(fill=NA),
    title=element_text(size=fsize,face="bold"),
    axis.title=element_text(size=fsize,face="bold"),
    axis.text = element_text(hjust=1,size = fsize, angle = 0),
    axis.line=element_line(colour="black"),
    panel.background=element_rect(fill = "white"),
    #    panel.background=element_blank(),
    #    panel.grid.major=element_line(linetype="dashed",colour="grey",size=0.5),
    panel.grid.major=element_blank(),
    strip.background=element_rect(fill="white", colour="white"),
    strip.text.x = element_text(size=fsize, colour = "black", angle = 0,face="bold"),
    legend.title = element_text(size=fsize),
    legend.text = element_text(size=fsize)
  )

#Defining function
#Making two figures
fourfigure <- function(p1,p2,p3,p4,x){
  grid.newpage()
  l <- grid.layout( nrow=2, ncol=2, heights=unit(c(1, 1.25),c("null","null")), widths=unit(c(1, 1),c("null","null")))
  png(filename=x, width=2400, height=1500,res=200)
  v <- viewport(layout = l) # [l]を作図領域として[v]に格納する
  pushViewport(v) #グラフを埋め込む前に，pushViewport()関数に作図領域を指定する。
  print(p1, vp = viewport(layout.pos.row=1,layout.pos.col=1))
  print(p2, vp = viewport(layout.pos.row=2,layout.pos.col=1))
  print(p3, vp = viewport(layout.pos.row=1,layout.pos.col=2))
  print(p4, vp = viewport(layout.pos.row=2,layout.pos.col=2))
  popViewport()  #作図を終了する
  dev.off()  
  return()
}
nc <- 0
xhist <- 0
Dataextract <- function(x,y,v,reg,cr,yr){
  xhist <- x[x$VAR==v & x$Year==yr & x$Region==reg & x$SSP=="SSP2" & x$Crop==cr, c(-1,-2,-3,-4,-10)]
  nocc1 <- y[y$VAR==v & y$Year==yr & y$Region==reg & y$SSP=="SSP2" & y$Crop==cr, c(6)]  
  names(xhist) <- c("GCM","RCP","CO2fert","SYR","N","Value")	#Rename the column
  minv <- min(min(xhist[6]),nocc1)
  if (varname=="PPOP_UNSH"){
  maxv <- max(max(xhist[6]),nocc1)
    bw <- maxv/100
  } else {
    maxv <- max(max(xhist[6]),nocc1)
    bw <- diff(range(xhist$Value))/100
  }
  
  mainlabel <- paste(reg,cr,yr)
  interact<-interaction(xhist$RCP, xhist$CO2fert, sep=" : ")
  
  mu <- ddply(xhist, "interact", summarise, grp.mean=mean(Value))
  v5 <- ddply(xhist,"interact",summarise, grp.quantile=quantile(Value,0.05))
  plotx <- ggplot() + 
    geom_density(data=xhist, aes(x=Value,fill=interact,y=..density..,color=interact),alpha = 0,size=2, position = "identity") +  
    geom_histogram(data=xhist, aes(x=Value,fill=interact,y=..density..,color=interact),alpha = 0.3, position = "identity") +
    #    geom_histogram(alpha = 0.3, position = "identity",  binwidth=bw, color="black") +
    xlim(minv-abs(minv)*0.01, maxv+abs(maxv)*0.01) +
    ggtitle(mainlabel) +
    ylab("density") + xlab(labelname) +
    MyThemeLine +
    geom_vline(data=mu,aes(xintercept=grp.mean, color=interact), linetype="dashed", size=1) +
  #  geom_vline(data=v5,aes(xintercept=grp.quantile, color=interact), linetype="solid", size=1) +
    #  geom_text(data=mu,aes(x=grp.mean,y=0.023,label=grp.mean),angle=-90,size=5,vjust=-bw/10,hjust=-0.0001) +
    # geom_text(data=v5,aes(x=grp.quantile,y=0.023,label=grp.quantile),angle=-90,size=5,vjust=-bw/10,hjust=-0.0001) +
  geom_vline(aes(xintercept=nocc1), color="black", size=1.2) 
  #  geom_text(aes(x=nocc1,y=0.023,label=nocc1),angle=-90,size=5,vjust=-bw/10,hjust=-0.0001)
  #  geom_area(mapping = aex(y = ifelse(y<0.2)))
#  ggsave(plot1,filename=filenamehist, dpi=250,width=8, height=6)
  return(plotx)  
}

twofigure <- function(p1,p2,x){
  grid.newpage()
  l <- grid.layout( nrow=1, ncol=2, heights=unit(c(1, 1.25),c("null","null")), widths=unit(c(1, 1),c("null","null")))
  png(filename=x, width=2400, height=1500,res=200)
  v <- viewport(layout = l) # [l]を作図領域として[v]に格納する
  pushViewport(v) #グラフを埋め込む前に，pushViewport()関数に作図領域を指定する。
  print(p1, vp = viewport(layout.pos.row=1,layout.pos.col=1))
  print(p2, vp = viewport(layout.pos.row=1,layout.pos.col=2))
  popViewport()  #作図を終了する
  dev.off()  
  return()
}

#Making five figures
fivefigure <- function(p1,p2,p3,p4,p5,x,l,wd,ht,rs,x2){
  grid.newpage()
  png(filename=x, width=wd, height=ht,res=rs)
  v <- viewport(layout = l)
  pushViewport(v)
  grid.text(x2, vp = viewport(layout.pos.row = 1, layout.pos.col = 1:2))
  print(p1, vp = viewport(layout.pos.row=2,layout.pos.col=1))
  print(p2, vp = viewport(layout.pos.row=2,layout.pos.col=2))
  print(p3, vp = viewport(layout.pos.row=3,layout.pos.col=1))
  print(p4, vp = viewport(layout.pos.row=3,layout.pos.col=2))
  print(p5, vp = viewport(layout.pos.row=4,layout.pos.col=1))
  #  print(title=x2)
  popViewport()
  dev.off()  
  return()
}
#---- making directory
dirname0 <- paste('../output/png/', sep = '')
dir.create(dirname0)

var.in <- c("PCAL","PPOP_UNSH","YEXO","XPRP")

for(k in 1:length(var.in)){
  varname <- var.in[k]
  
  #  dirname <- paste('../output/png/', varname, sep = '')
  dirname <- paste(dirname0, varname, sep = '')
  dir.create(dirname)
  
}
#--- data load
yearname <- c('2050')
pdata <- rgdx.param('../data/Dataall_202050.gdx','Data')
pdata <- rgdx.param('../data/Dataall_202050.gdx','Data_agg')
names(pdata) <- colnames
NoCC_load <- rgdx.param('../data/Dataall_202050.gdx','Data_Nocc')
names(NoCC_load) <- colnamesNocc

pdata$RCP <- factor(pdata$RCP, levels=c("RCP85","RCP26"))
pdata$Region <- factor(pdata$Region, levels=c("USA","CAN","JPN","TUR","XER","XE25","XOC","CIS","XNF","CHN","BRA","IND","XLM","XSE","XME","XSA","XAF","WLD"))
pdata$Region <- factor(pdata$Region, levels=c("Other","US","EU","FSU","CHN","BRA","China","LAM","SEAsia","SAsia","Afr_ME","WLD"))

fao <- rgdx.param('../data/Risk_fao.gdx','risk_fao')
fao <- rgdx.param('../data/Risk_fao.gdx','pcal_fao')
names(fao) <- c("Reg","Year","Value")
faoworld <- fao %>% dplyr::filter(Reg == "WLD")


#Figure 1-c, Figure2
#-----#1 PPOP_UNSH & PCAL: CDF Histogram for a single region ----
var.in <- c("PCAL","PPOP_UNSH")
var.in <- c("PPOP_UNSH")
var.in <- c("YEXO")
var.in <- c("PCAL","PPOP_UNSH","XPRP","YEXO")
var.in <- c("XPRP","YEXO")
regionlist <- c("IND","XSA","XAF","BRA","CIS","XNF","CHN","XLM","XSE","XME")
regionlist <- c("USA","CAN","JPN","TUR","XER","XE25","XOC","CIS","XNF","CHN","BRA","IND","XLM","XSE","XME","XSA","XAF","WLD")
regionlist <- c("XSA","XAF","WLD")
regionlist <- c("WLD")
regionlist <- c("IND")

crop_list <- c("PDR","WHT","GRO","OSD") 
crop_list <- c("PDR") 
nocc <- 0
cropname <- c("dum")
yearlist <- c("2050")
ssp_list <- c("SSP2","SSP3")
ssp_list <- c("SSP2")
ssp_list <- c("SSP3")
x <- pdata
y <- NoCC_load

yrn <- 1
k <- 1
r <- 1


for(yrn in 1:length(yearlist)){
  yearname <- yearlist[yrn]
  for(k in 1:length(var.in)){
    varname <- var.in[k]
    
    for(r in 1:length(regionlist)){
      regionname <- regionlist[r]
      
      if (varname=="PCAL"){
        labelname <- c("Calorie intake [kcal/person/day]")
        cropnum <- 1
      } else if (varname=="PPOP_UNSH"){
        labelname <- c("Population at risk of hunger [billion]")
        cropnum <- 1
      } else if (varname=="YEXO"){
        labelname <- c("Yield [tonne/ha]")
        cropnum <- 4
      } else if (varname=="XPRP"){
        labelname <- c("Price index [-]")
        cropnum <- 4
      } else {
        labelname <- c("aaaaaaaa")
      }
      
            if (cropnum==4){
#        crop_list <- c("PDR","WHT","GRO","OSD")
         crop_list <- c("WHT","GRO","OSD","PDR") 
      }else if (cropnum==1){   
        crop_list <- c("dum")
      }
      for (crn in 1:length(crop_list)){
        cropname <- crop_list[crn]
      
        if (varname=="PCAL"){
          mainlabel <- paste(regionname)
        } else if (varname=="PPOP_UNSH"){
          mainlabel <- paste(regionname)
        } else if (varname=="YEXO"){
          mainlabel <- paste(regionname,cropname)
        } else if (varname=="XPRP"){
          mainlabel <- paste(regionname,cropname)
        } else {
          mainlabel <- c("aaaaaaaa")
        }
          
        for (ssp in 1:length(ssp_list)){
          flag <- 0
          if (ssp_list[ssp]=="SSP2"){
            xhist <- x[x$VAR==varname & x$Year==yearname & x$Region==regionname & x$SSP=="SSP2" & x$Crop==cropname, c(-1,-2,-3,-4,-10)]
            nocc1 <- y[y$VAR==varname & y$Year==yearname & y$Region==regionname & y$SSP=="SSP2" & y$Crop==cropname, c(6)]  
            names(xhist) <- c("GCM","RCP","XXX","SYR","N","Value")	#Rename the column
            flag <-1
          }else if (ssp_list[ssp]=="SSP3"){
            xhist <- x[x$VAR==varname & x$Year==yearname & x$Region==regionname & x$SSP=="SSP3" & x$Crop==cropname, c(-1,-2,-3,-4,-10)]
            nocc1 <- y[y$VAR==varname & y$Year==yearname & y$Region==regionname & y$SSP=="SSP3" & y$Crop==cropname, c(6)]  
            names(xhist) <- c("GCM","RCP","XXX","SYR","N","Value")	#Rename the column
            flag <-1
          }else{
            if (yearname=="2050"){
              xhist <- x[x$VAR==varname & x$Year==yearname & x$Region==regionname & x$CO2fert=="noco2" & x$Crop==cropname, c(-1,-2,-3,-7,-10)]
              nocc1 <- y[y$VAR==varname & y$Year==yearname & y$Region==regionname & y$Crop==cropname, c(4,6)]  
              names(xhist) <- c("XXX","GCM","RCP","SYR","N","Value")	#Rename the column
              flag <-1
            }
          }
          if (flag != 0){
            
            minv <- min(min(xhist[6]),min(nocc1[2]))
            minv <- min(min(xhist[6]),nocc1)
            
            if (varname=="PPOP_UNSH"){
              maxv <- max(max(xhist[6]),min(nocc1[2]))
              maxv <- max(max(xhist[6]),nocc1)
              bw <- maxv/100
            } else {
              maxv <- max(max(xhist[6]),min(nocc1[2]))
              bw <- diff(range(xhist$Value))/100
            }


            
            interact<-interaction(xhist$RCP, xhist$XXX, sep=" : ")
            xhist <- cbind(interact,xhist)
              
            col<- c('mediumblue','orangered3','mediumaquamarine','orange')
            col<- c('blue','red','cyan','orange')
          mu <- ddply(xhist,"interact", summarise, grp.mean=mean(Value))
          v5 <- ddply(xhist,"interact",summarise, grp.quantile=quantile(Value,0.05))
          labelname <- paste(labelname)
          labelnamelog <- paste('log10 (',labelname,')')
        
          plotx <- ggplot() + 
            geom_step(data=xhist, aes(x=Value,y=..y..,color=interact),stat="ecdf") +
            scale_x_log10() +
#            xlim(minv-abs(minv)*0.01, maxv+abs(maxv)*0.01) +
            scale_fill_manual(name="Impact range",values=col) +
            scale_color_manual(name="Impact range",values=col) +
            ggtitle(mainlabel) +
#            ggtitle("SSP3") + theme(plot.title=element_text(hjust=1.5)) +
            ylab("Density") + 
            xlab(labelnamelog) +
            MyThemeLine +
            geom_vline(data=mu,aes(xintercept=grp.mean, color=interact), linetype="dashed", size=1) +
            #  geom_vline(data=v5,aes(xintercept=grp.quantile, color=interact), linetype="solid", size=1) +
            #  geom_text(data=mu,aes(x=grp.mean,y=0.023,label=grp.mean),angle=-90,size=5,vjust=-bw/10,hjust=-0.0001) +
            # geom_text(data=v5,aes(x=grp.quantile,y=0.023,label=grp.quantile),angle=-90,size=5,vjust=-bw/10,hjust=-0.0001) +
#            geom_vline(aes(xintercept=nocc1), color="black", size=1.2) 
#            theme(legend.position="none")
          geom_hline(yintercept=0.1, color="black", linetype = "dashed") +
          geom_hline(yintercept=0.01, color="black", linetype = "dashed")
      
            
          plot(plotx)
          filenamehist <- paste('../output/png/CDF/', varname ,'_',ssp_list[ssp],regionname,'_',yearname,cropname,'.png', sep = '')
          
          ggsave(plotx,filename=filenamehist, dpi=250,width=8.5, height=6)
          }
        }
      }
    }
  }
}  
#----End --- PPOP_UNSH & PCAL: Histogram for a single region----

#-----#1 PPOP_UNSH & PCAL: CDF by GCMs for a single region ----
var.in <- c("PCAL","PPOP_UNSH")
var.in <- c("PPOP_UNSH")
var.in <- c("YEXO")
var.in <- c("PCAL","PPOP_UNSH","XPRP","YEXO")
var.in <- c("XPRP","YEXO")
regionlist <- c("IND","XSA","XAF","BRA","CIS","XNF","CHN","XLM","XSE","XME")
regionlist <- c("USA","CAN","JPN","TUR","XER","XE25","XOC","CIS","XNF","CHN","BRA","IND","XLM","XSE","XME","XSA","XAF","WLD")
regionlist <- c("XSA","XAF","WLD")
regionlist <- c("WLD")
regionlist <- c("IND")

crop_list <- c("PDR","WHT","GRO","OSD") 
crop_list <- c("PDR") 
nocc <- 0
cropname <- c("dum")
yearlist <- c("2050")
ssp_list <- c("SSP2","SSP3")
ssp_list <- c("SSP2")
ssp_list <- c("SSP3")
x <- pdata
y <- NoCC_load

yrn <- 1
k <- 1
r <- 1

gcm.names <- c("mes_","ipsl","ne1m","ge2m","hg2e")
gcm.in <-  gcm.names


for(yrn in 1:length(yearlist)){
  yearname <- yearlist[yrn]
  for(k in 1:length(var.in)){
    varname <- var.in[k]
    
    for(r in 1:length(regionlist)){
      regionname <- regionlist[r]
      
      if (varname=="PCAL"){
        labelname <- c("Calorie intake [kcal/person/day]")
        cropnum <- 1
      } else if (varname=="PPOP_UNSH"){
        labelname <- c("Population at risk of hunger [billion]")
        cropnum <- 1
      } else if (varname=="YEXO"){
        labelname <- c("Yield [tonne/ha]")
        cropnum <- 4
      } else if (varname=="XPRP"){
        labelname <- c("Price index [-]")
        cropnum <- 4
      } else {
        labelname <- c("aaaaaaaa")
      }
      
      if (cropnum==4){
        #        crop_list <- c("PDR","WHT","GRO","OSD")
        crop_list <- c("WHT","GRO","OSD","PDR") 
      }else if (cropnum==1){   
        crop_list <- c("dum")
      }
      for (crn in 1:length(crop_list)){
        cropname <- crop_list[crn]
        
        
        for (ssp in 1:length(ssp_list)){
          flag <- 0
          if (ssp_list[ssp]=="SSP2"){
            xhist <- x[x$VAR==varname & x$Year==yearname & x$Region==regionname & x$SSP=="SSP2" & x$Crop==cropname, c(-1,-2,-3,-4,-10)]
            nocc1 <- y[y$VAR==varname & y$Year==yearname & y$Region==regionname & y$SSP=="SSP2" & y$Crop==cropname, c(6)]  
            names(xhist) <- c("GCM","RCP","XXX","SYR","N","Value")	#Rename the column
            flag <-1
          }else if (ssp_list[ssp]=="SSP3"){
            xhist <- x[x$VAR==varname & x$Year==yearname & x$Region==regionname & x$SSP=="SSP3" & x$Crop==cropname, c(-1,-2,-3,-4,-10)]
            nocc1 <- y[y$VAR==varname & y$Year==yearname & y$Region==regionname & y$SSP=="SSP3" & y$Crop==cropname, c(6)]  
            names(xhist) <- c("GCM","RCP","XXX","SYR","N","Value")	#Rename the column
            flag <-1
          } else {
            if (yearname=="2050"){
              xhist <- x[x$VAR==varname & x$Year==yearname & x$Region==regionname & x$CO2fert=="noco2" & x$Crop==cropname, c(-1,-2,-3,-7,-10)]
              nocc1 <- y[y$VAR==varname & y$Year==yearname & y$Region==regionname & y$Crop==cropname, c(4,6)]  
              names(xhist) <- c("XXX","GCM","RCP","SYR","N","Value")	#Rename the column
              flag <-1
            }
          }
          if (flag != 0){
            
            minv <- min(min(xhist[6]),min(nocc1[2]))
            minv <- min(min(xhist[6]),nocc1)
            
            if (varname=="PPOP_UNSH"){
              maxv <- max(max(xhist[6]),min(nocc1[2]))
              maxv <- max(max(xhist[6]),nocc1)
              bw <- maxv/100
            } else {
              maxv <- max(max(xhist[6]),min(nocc1[2]))
              bw <- diff(range(xhist$Value))/100
            }
            
            
            
            interact<-interaction(xhist$RCP, xhist$XXX, sep=" : ")
            xhist <- cbind(interact,xhist)
            
            col<- c('mediumblue','orangered3','mediumaquamarine','orange')
            col<- c('blue','red','cyan','orange')
            mu <- ddply(xhist,"interact", summarise, grp.mean=mean(Value))
            v5 <- ddply(xhist,"interact",summarise, grp.quantile=quantile(Value,0.05))
            labelname <- paste(labelname)
            labelnamelog <- paste('log10 (',labelname,')')
            
            for(gcm in 1:length(gcm.in)){
              #gcm <- 1
              gcmname <- gcm.in[gcm]
  
              xhist2 <- xhist[xhist$GCM==gcmname, c(-2)]  
              
              if (varname=="PCAL"){
                mainlabel <- paste(regionname)
              } else if (varname=="PPOP_UNSH"){
                mainlabel <- paste(regionname)
              } else if (varname=="YEXO"){
                mainlabel <- paste(regionname,cropname,gcmname)
              } else if (varname=="XPRP"){
                mainlabel <- paste(regionname,cropname)
              } else {
                mainlabel <- c("aaaaaaaa")
              }
              
              
            plotx <- ggplot() + 
              geom_step(data=xhist2, aes(x=Value,y=..y..,color=interact),stat="ecdf") +
              scale_x_log10() +
              #            xlim(minv-abs(minv)*0.01, maxv+abs(maxv)*0.01) +
              scale_fill_manual(name="Impact range",values=col) +
              scale_color_manual(name="Impact range",values=col) +
              ggtitle(mainlabel) +
              #            ggtitle("SSP3") + theme(plot.title=element_text(hjust=1.5)) +
              ylab("Density") + 
              xlab(labelnamelog) +
              MyThemeLine +
              geom_vline(data=mu,aes(xintercept=grp.mean, color=interact), linetype="dashed", size=1) +
              #  geom_vline(data=v5,aes(xintercept=grp.quantile, color=interact), linetype="solid", size=1) +
              #  geom_text(data=mu,aes(x=grp.mean,y=0.023,label=grp.mean),angle=-90,size=5,vjust=-bw/10,hjust=-0.0001) +
              # geom_text(data=v5,aes(x=grp.quantile,y=0.023,label=grp.quantile),angle=-90,size=5,vjust=-bw/10,hjust=-0.0001) +
              #            geom_vline(aes(xintercept=nocc1), color="black", size=1.2) 
              #            theme(legend.position="none")
              geom_hline(yintercept=0.1, color="black", linetype = "dashed") +
              geom_hline(yintercept=0.01, color="black", linetype = "dashed")
            
            
            plot(plotx)
            filenamehist <- paste('../output/png/CDF/', varname ,'_',ssp_list[ssp],regionname,'_',yearname,cropname,'_gcm.png', sep = '')
          # ggsave(plotx,filename=filenamehist, dpi=250,width=8.5, height=6)
            
            if (gcm==1){
              g1 <- plotx
            }else if(gcm==2){
              g2 <- plotx
            }else if(gcm==3){
              g3 <- plotx
            }else if(gcm==4){
              g4 <- plotx
            }else if(gcm==5){
              g5 <- plotx
            }
            
            }  
            ht1 <- 3000
            rs1 <- 500
            wd1 <- 5000
            l <- grid.layout(nrow=4, ncol=2, widths=unit(c(1.5, 1.5),c("null","null")), heights=unit(c(0.1, 1, 1,1),c("null","null","null","null")))
            file_name_out <- filenamehist
            title <- paste('YEXO',regionname,cropname)
            fivefigure(g1,g2,g3,g4,g5,file_name_out,l,wd1,ht1,rs1,title)
            
            plot(g5)
            
            dev.off()
            
          }
        }
      }
    }
  }
}  
#----End --- PPOP_UNSH & PCAL: CDF by GCMs for a single region----
# Figure1-a
#-----#2 PPOP_UNSH & PCAL: Ribbon for a single region ----
var.in <- c("PCAL")
var.in <- c("PPOP_UNSH")
regionlist <- c("WLD","IND","XSA","XAF","BRA")
regionlist <- c("WLD")

crop_list <- c("PDR","WHT","GRO","OSD") 
nocc <- 0
cropname <- c("dum")
yearlist <- c("2020","2050")
ssp_list <- c("SSP2","SSP3")
x <- pdata
y <- NoCC_load


#######
####### loading long term analysis results
#######

pdata_long <- rgdx.param('../data/Dataall_longanalysis.gdx','Pagmipout')

nocct <- pdata_long[pdata_long$R=="WLD" & pdata_long$Sgcm=="nocc" & pdata_long$Srcp=="nocc" &  
                      pdata_long$co2fert=="nocc" & ( pdata_long$Title_Agmip=="CALO" | pdata_long$Title_Agmip=="PPOP_UNSH"), c(1,5,8,10)]
colnamesnocct <- c("SSP","Year","VAR","Value")
names(nocct) <- colnamesnocct
levels(nocct$VAR)[levels(nocct$VAR)=="CALO"] <- "PCAL"
nocct[2] <- lapply(nocct[2],as.integer)
nocct[2] <- (nocct[2]-1)*5+2005
nocct <- nocct[(nocct$Year>=2020 & nocct$Year<=2050), c(1,2,3,4)]

cct <- pdata_long[pdata_long$R=="WLD" & pdata_long$Srcp=="RCP85" &  
                      pdata_long$co2fert=="noco2" & ( pdata_long$Title_Agmip=="CALO" | pdata_long$Title_Agmip=="PPOP_UNSH"), c(1,2,5,8,10)]
colnamescct <- c("SSP","GCM","Year","VAR","Value")
names(cct) <- colnamescct
levels(cct$VAR)[levels(cct$VAR)=="CALO"] <- "PCAL"
cct[3] <- lapply(cct[3],as.integer)
cct[3] <- (cct[3]-1)*5+2005
cct <- cct[(cct$Year>=2020 & cct$Year<=2050), c(1,2,3,4,5)]

#######
#######
#######


for(k in 1:length(var.in)){
  varname <- var.in[k]
  if (varname=="PCAL"){
    labelname <- c("Calorie intake [kcal/person/day]")
    cropnum <- 1
  } else if (varname=="PPOP_UNSH"){
    labelname <- c("Population at risk of hunger [billion]")
    cropnum <- 1
  } else if (varname=="YEXO"){
    labelname <- c("Yield [tonne/ha]")
    cropnum <- 4
  } else if (varname=="XPRP"){
    labelname <- c("Price index [-]")
    cropnum <- 4
  } else {
    labelname <- c("aaaaaaaa")
  }
  for(r in 1:length(regionlist)){
    regionname <- regionlist[r]
    
    if (cropnum==4){
      crop_list <- c("PDR","WHT","GRO","OSD") 
    }else if (cropnum==1){   
      crop_list <- c("dum")
    }
    for (crn in 1:length(crop_list)){
      cropname <- crop_list[crn]
      #        for (ssp in 1:length(ssp_list)){
      flag <- 0
      xrib <- x[x$VAR==varname & x$RCP=="RCP85" & x$Region==regionname & x$Crop==cropname & x$CO2fert=="noco2", c(-1,-3,-6,-10)]
      #            nocc1 <- y[y$VAR==varname & y$Region==regionname & y$SSP=="SSP2" & y$Crop==cropname, c(6)]  
      names(xrib) <- c("Year","SSP","GCM","XXX","SYR","N","Value")	#Rename the column
      
      minv <- min(min(xrib[7]))
      if (varname=="PPOP_UNSH"){
        maxv <- max(max(xrib[7]))
        bw <- maxv/100
      } else {
        maxv <- max(max(xrib[7]))
        bw <- diff(range(xrib$Value))/100
      }
      
      mainlabel <- paste(regionname,cropname)
      faoworld2 <- faoworld
      faoworld2[2] <- lapply(faoworld2[2],as.integer)
      faoworld2[2] <- faoworld2[2]+1990
      faoworld2 <- faoworld2[(faoworld2$Year<=2005), c(1,2,3)]
      faoworld2 <- faoworld2[(faoworld2$Year<=2010), c(1,2,3)]
      
      #      interact<-interaction(xrib$RCP, xrib$SSP, sep=" : ")
      
      xrib1 <- ddply(xrib,c("Year","SSP"),transform,max=max(Value))
      xrib2 <- ddply(xrib1,c("Year","SSP"),transform,min=min(Value))
      xrib5 <- ddply(xrib2,c("Year","SSP"),transform,v1=quantile(Value,0.01))
      xrib6 <- ddply(xrib5,c("Year","SSP"),transform,v2=quantile(Value,0.02))
      xrib7 <- ddply(xrib6,c("Year","SSP"),transform,v5=quantile(Value,0.05))
      xrib8 <- ddply(xrib7,c("Year","SSP"),transform,v10=quantile(Value,0.1))
      xrib9 <- ddply(xrib8,c("Year","SSP"),transform,v20=quantile(Value,0.2))
      xrib10 <- ddply(xrib9,c("Year","SSP"),transform,v33=quantile(Value,0.3333))
      xrib105 <- ddply(xrib10,c("Year","SSP"),transform,v50=quantile(Value,0.5))
      xrib11 <- ddply(xrib105,c("Year","SSP"),transform,v66=quantile(Value,0.6666))
      xrib12 <- ddply(xrib11,c("Year","SSP"),transform,v80=quantile(Value,0.8))
      xrib13 <- ddply(xrib12,c("Year","SSP"),transform,v90=quantile(Value,0.9))
      xrib14 <- ddply(xrib13,c("Year","SSP"),transform,v95=quantile(Value,0.95))
      xrib15 <- ddply(xrib14,c("Year","SSP"),transform,v98=quantile(Value,0.98))
      xrib16 <- ddply(xrib15,c("Year","SSP"),transform,v99=quantile(Value,0.99))
      
      names(xrib16)[8]<-paste("max")
      names(xrib16)[9]<-paste("min")
      
      #            col<- c('#00FF00', '#0000FF',"#000000")
      col<- c('mediumblue','orangered3',"#000000")
      
      xrib16[1] <- lapply(xrib16[1],as.integer)
      xrib16[1] <- (xrib16[1]-18)*10+2010
      xrib17 <- xrib16 %>% dplyr::distinct(Year,SSP,GCM,XXX,max,min,v1,v2,v5,v10,v20,v33,v50,v66,v80,v90,v95,v98,v99)
#      xrib17 <- xrib16 %>% 
#        distinct(Year,SSP,GCM,XXX,max,min,v1,v2,v5,v10,v20,v33,v50,v66,v80,v90,v95,v98,v99)
      #2015 continuation
      fao2015 <- dplyr::filter(faoworld,Year==2005)
      fao2015 <- dplyr::filter(faoworld,Year==2010)

      fao2015 <- as.numeric(fao2015[1,3])
      xrib16add <- xrib17 %>%
        dplyr::mutate(v000=fao2015,v00=fao2015,v01=fao2015,v1=fao2015,v0=fao2015,v2=fao2015,v3=fao2015,v4=fao2015,v5=fao2015,v6=fao2015,v7=fao2015,v8=fao2015,v9=fao2015,v10=fao2015,v11=fao2015,v12=fao2015,v13=fao2015,v14=fao2015,v15=fao2015) %>%
        dplyr::select(-max,-min,-v1,-v2,-v5,-v10,-v20,-v33,-v50,-v66,-v80,-v90,-v95,-v98,-v99)
      names(xrib16add) <- c("Year","SSP","GCM","XXX","max","min","v1","v2","v5","v10","v20","v33","v50","v66","v80","v90","v95","v98","v99")
      xrib16add[1] <- 2005
      xrib16add[1] <- 2010
      xrib18 <- rbind(xrib17,xrib16add)

      cct0 <- cct[cct$VAR==varname,c(1,2,3,4,5)]
      cct1 <- ddply(cct0,c("Year","SSP"),transform,max=max(Value))
      cct2 <- ddply(cct1,c("Year","SSP"),transform,min=min(Value))
  
      nocct_risk <- nocct[nocct$VAR==varname, c(1,2,4)]

      nocct_base <- fao2015
      nocct_base[2] <- 2005
      nocct_base[2] <- 2010
      nocct_base_SSP2 <- nocct_base
      nocct_base_SSP2[3] <- c("SSP2")
      nocct_base_SSP2 <- nocct_base_SSP2[c(3,2,1)]
      nocct_base_SSP3 <- nocct_base
      nocct_base_SSP3[3] <- c("SSP3")
      nocct_base_SSP3 <- nocct_base_SSP3[c(3,2,1)]
      nocct_base_SSP2[2] <- 2010
      nocct_base_SSP3[2] <- 2010    
      nocct_base_SSP2[2] <- 2005
      nocct_base_SSP3[2] <- 2005    
      nocct_risk1<- rbind(nocct_base_SSP2,nocct_base_SSP3,nocct_risk)
      nocct_risk1[2] <- lapply(nocct_risk1[2],as.integer)
      nocct_risk1[3] <- lapply(nocct_risk1[3],as.numeric)
      
            plotrib <-ggplot()+ 
        #        geom_ribbon(data=xrib18,aes(x=Year,ymin=min, ymax=max, fill = SSP), alpha=0.1) +
        geom_ribbon(data=xrib18,aes(x=Year,ymin=v1, ymax=v99, fill = SSP), alpha=0.2) +
        geom_ribbon(data=xrib18,aes(x=Year,ymin=v5, ymax=v95, fill = SSP), alpha=0.3) +
        geom_ribbon(data=xrib18,aes(x=Year,ymin=v10, ymax=v90, fill = SSP), alpha=0.4) +
        geom_ribbon(data=xrib18,aes(x=Year,ymin=v33, ymax=v66, fill = SSP), alpha=0.7) +
        geom_line(data=xrib18,aes(x=Year,y=v50,group = SSP, colour=SSP), size = 2) +
        labs(aes(colour = SSP)) + 
        scale_colour_manual(name="Median value",values=col) + 
        scale_fill_manual(name="Impact range",values=col) +
        theme_bw(base_size = 30) + labs(y =labelname, x="") +
        #              ylim(minv-abs(minv)*0.01, maxv+abs(maxv)*0.01) +
        MyThemeLine +
        #        ggtitle(mainlabel) +
        geom_line(data=faoworld2,aes(x=Year,y=Value,group=Reg)) +
        geom_line(data=nocct_risk1,aes(x=Year,y=Value,fill=SSP),size=1.0) +
        geom_line(data=cct2,aes(x=Year,y=max,fill=SSP),size=1.0,linetype = 2) +
        geom_line(data=cct2,aes(x=Year,y=min,fill=SSP),size=1.0,linetype = 2)  
#        ylim(2500,3000)             
      #        xlim(1990,2050)
      #        ylim(minv-abs(minv)*0.01, maxv+abs(maxv)*0.01) 
      
      plot(plotrib)
      
      filenamerib <- paste('../output/png/', varname ,'/',regionname,'_',yearname,cropname,'ribbon_rcp85v2.png', sep = '')
      ggsave(plotrib,filename=filenamerib, dpi=250,width=8, height=6)
      #        }
    }
  }
}
#Figure3, 1-b
#----#3 return period
ssp_list <- c("SSP2")
rcp_list <- c("RCP85")
var.in <- c("PPOP_UNSH")
allregion <- c("WLD")
co2_list <- c("noco2")
yearlist <- c("2050")
regionlist <- c("WLD")


ssp_list <- c("SSP2","SSP3")
var.in <- c("PCAL")
var.in <- c("PCAL","PPOP_UNSH")
varname <- var.in
rcp_list <- c("RCP85","RCP26")
allregion <- c("WLD","USA","XE25","XER","TUR","XOC","CHN","IND","JPN","XSE","XSA","CAN","BRA","XLM","CIS","XME","XNF","XAF")
co2_list <- c("co2","noco2")
co2_list <- c("noco2")
yearlist <- c("2020","2030","2040","2050")
regionlist <- c("WLD")


z <- pdata[pdata$Crop=="dum", c(-5,-8,-9,-10)]

#names(z) <- c("VAR","Year","Region","SSP","RCP","CO2fert","Value")	#Rename the column

yall <- 0
yall <- data.frame(1,1,1,1,1,1,1,1)	#Rename the column
names(yall) <- c("VAR","Year","Region","RCP","SSP","CO2fert","V","return")	#Rename the column

returnp<-c(0,0.01,0.02,0.05,0.1,0.2,0.33,0.5,0.66,0.8,0.9,0.95,0.98,0.99,1)
returnn<-c("MIN","RP001","RP002","RP005","RP01","RP02","RP033","RP05","RP066","RP08","RP09","RP095","RP098","RP099","MAX")

for(k in 1:length(var.in)){
for(yr in 1:length(yearlist)){
 for (ssp in 1 :length(ssp_list)){
   for (co2 in 1 : length(co2_list)){
     for (rc in 1 : length(rcp_list)){
        for(r in 1:length(allregion)){
          for(i in 1 : 15){

                        fl <- z[z$VAR==var.in[k] & z$Year==yearlist[yr] & z$SSP==ssp_list[ssp] & z$CO2fert==co2_list[co2] & z$RCP==rcp_list[rc] & z$Region==allregion[r],c(7)]
            if (length(fl)>0){
              a <- returnp[i]
              b <- returnn[i]
              v1 <- quantile(fl,a)
              y5 <- data.frame(var.in[k],yearlist[yr],allregion[r],rcp_list[rc],ssp_list[ssp],co2_list[co2],b,v1[1])	#Rename the column
              names(y5) <- c("VAR","Year","Region","RCP","SSP","CO2fert","V","return")	#Rename the column
                yall <- rbind(yall,y5) 
            y5<-0
            }
          }
        }
      }
    }
  }
}
}
symDim <- 8
attr(yall, "symName") <- "Foodshock"
lst <- wgdx.reshape(yall, symDim)
wgdx.reshape(yall, symDim, gdxName = "../data/FoodshockR.gdx")
system(paste("gams","../prog/foodstock.gms",sep=" "))

ReturnPeriodFS <- rgdx.param('../output/gdx/PVAR.gdx','PVARagg')
ReturnPeriodpr <- rgdx.param('../output/gdx/PVAR.gdx','PVARRPagg')
ReturnPeriodpr <- ReturnPeriodpr[c(8)]
#ReturnPeriodFS[7] <- ReturnPeriodFS[7]*100
ReturnPeriod0 <- cbind(ReturnPeriodFS,ReturnPeriodpr)
colnamesRP0 <- c("VAR","Year","Region","RCP","SSP","CO2fert","RetPerI","Stock","RetPer")
names(ReturnPeriod0) <- colnamesRP0

ReturnPeriodSDM <- rgdx.param('../output/gdx/PVAR.gdx','PSDMagg')
ReturnPeriodSDMpr <- rgdx.param('../output/gdx/PVAR.gdx','PSDMRPagg')
ReturnPeriodSDMpr <- ReturnPeriodSDMpr[c(8)]
SDM0 <- cbind(ReturnPeriodSDM,ReturnPeriodSDMpr)
names(SDM0) <- colnamesRP0

actualstock <- rgdx.param('../output/gdx/PVAR.gdx','ActualStock')
colnamesAS <- c("Region","Stock")
names(actualstock) <- colnamesAS


for(k in 1:length(var.in)){
  varname <- var.in[k]
  
  colnamesRP <- c("Year","Region","RCP","SSP","CO2fert","RetPerI","Stock","RetPer")
  ReturnPeriod  <- ReturnPeriod0[ReturnPeriod0$VAR==varname, c(2,3,4,5,6,7,8,9)]
  names(ReturnPeriod) <- colnamesRP
 
  SDM  <- SDM0[SDM0$VAR==varname, c(2,3,4,5,6,7,8,9)]
  names(SDM) <- colnamesRP
  
  if (varname=="PCAL"){
    ylabelname <- c("Required food stock (EJ)")
  } else if (varname=="PPOP_UNSH"){
    ylabelname <- c("Change in population at risk of hunger\n [billion]")
  } else {
    ylabelname <- c("aaaaaaaa")
  }
  

for(yr in 1:length(yearlist)){
  yearname <-yearlist[yr]
  for (ssp in 1 :length(ssp_list)){
    sspname <- ssp_list[ssp]
    
x0 <- ReturnPeriod[ReturnPeriod$Year==yearname & ReturnPeriod$CO2fert=="noco2" & ReturnPeriod$SSP==sspname,c(-1,-4,-5,-6)]
names(x0) <- c("Region","RCP","Stock","RetPer")
x <- regionlabelagg(x0)

plot1 <- ggplot() + 
  geom_point(data=x, aes(x=RetPer,y=Stock,group=RCP,color=RCP)) + MyThemeLine +
  xlab("return period (Year)")+ylab("Required food stock (EJ)")
plot2 <- plot1+ facet_wrap( ~ Region,scales="free")
filename <- paste("../output/png/return/",yearlist[yr],"return_stock.png")
ggsave(plot2,filename=filename, dpi=250,width=10, height=6)

plot3 <- ggplot() + 
  geom_area(data=x[x$Region!="WLD" & x$RCP=="RCP85",c(1,2,3,4)], aes(x=RetPer,y=Stock,fill = Region),stat = "identity") + MyThemeLine + scale_fill_manual(values =linepalette)+
  xlab("Return period (Year)")+ylab("Required food stock (EJ)") +
geom_line(data=x[x$Region=="WLD" & x$RCP=="RCP85",c(1,2,3,4)], aes(x=RetPer,y=Stock,color=Region),stat = "identity") + MyThemeLine + scale_color_manual(values =c("black")) +
  annotate("segment",x=0,xend=100,y=0,yend=0,linetype="solid",color="grey")

filename <- paste("../output/png/return/",yearlist[yr],sspname,"return_stock_reg.png")
plot(plot3)
ggsave(plot3,filename=filename, dpi=250,width=10, height=6)


# bar chart comparing actal & required stocks
reqstock <- x0[x0$Region!="WLD" & x0$RCP=="RCP85" & x0$RetPer==100,c(1,3)]
reqstock[3] <- c("Required \n stock")
actualstock1 <- actualstock[actualstock$Region!="WLD",c(1,2)]
actualstock1[3] <- c("Actual \n stock")

stocktable <- rbind(reqstock,actualstock1)
stocktable1 <- regionlabelagg(stocktable) 

gbar_stock <- ggplot() + 
  geom_bar(data=stocktable1, aes(x=as.factor(V3),y=Stock,fill=Region), position="stack", stat="identity", colour='black') +
  ylab("Food stock (EJ)") + xlab("") +	MyThemeLine +
  scale_fill_manual(values =linepalette) + theme(axis.text.x=element_text(hjust=0.5))
  
plot(gbar_stock)  
#filenamestack <- paste("../output/png/return/",yearname,sspname,varname,"return_regstack.png")
filenamestack <- paste("../output/png/return/",yearname,sspname,"return_regstack.png")
ggsave(gbar_stock,filename=filenamestack, dpi=250,width=7, height=6)

  }
}
x <- ReturnPeriod[ReturnPeriod$Year=="2050" & ReturnPeriod$CO2fert=="noco2" &  ReturnPeriod$RCP=="RCP85",c(-1,-3,-5,-6)]
names(x) <- c("Region","SSP","Stock","RetPer")
plot1 <- ggplot() + 
  geom_point(data=x, aes(x=RetPer,y=Stock,group=SSP,color=SSP)) + MyThemeLine +
  xlab("Return period (Year)")+ylab("Required food stock(EJ)")
plot2 <- plot1+ facet_wrap( ~ Region,scales="free") + ggtitle("RCP8.5")
filename <- paste("../output/png/return/",yearlist[yr],varname,"return_RCP85.png")
ggsave(plot2,filename=filename, dpi=250,width=10, height=6)

for (r in 1 :length(regionlist)){
 regionname <- regionlist[r]

#Figure1-b/c
x <- ReturnPeriod[(ReturnPeriod$Year=="2050" | ReturnPeriod$Year=="2020") & ReturnPeriod$Region==regionname &  ReturnPeriod$CO2fert=="noco2" &  
                    ReturnPeriod$RCP=="RCP85",c(-2,-3,-5,-6)]
names(x) <- c("Year","SSP","Stock","RetPer")
interact<-interaction(x$Year, x$SSP, sep=" : ")
xrp <- cbind(interact,x)
colrp <- c('royalblue','mediumblue','tomato','orangered3')
plot1 <- ggplot() + 
  geom_line(data=xrp, aes(x=RetPer,y=Stock,group=interact,color=interact), size=1.5) + MyThemeLine +
  scale_color_manual(values =colrp) +
  xlab("Return period (Year)")+
  ylab(ylabelname) +
  #+ ggtitle("RCP8.5") +
  geom_point(data=xrp, aes(x=RetPer,y=Stock,group=Year,shape=Year), size=3) +
  scale_shape_manual(values = c(15,1))
plot(plot1)
filename <- paste("../output/png/return/",yearname,regionname,varname,"return_SSP.png")
ggsave(plot1,filename=filename, dpi=250,width=8.5, height=6)

#Figure1-c/b SDM
xsdm <- SDM[ SDM$RetPer=="100" & SDM$Region==regionname &  SDM$CO2fert=="noco2" &  
           SDM$RCP=="RCP85",c(-2,-3,-5,-6,-8)]
names(xsdm) <- c("Year","SSP","Stock")
plot2 <- ggplot() + 
  geom_line(data=xsdm, aes(x=Year,y=Stock,group=SSP,color=SSP), size=1.5) + MyThemeLine +
  scale_color_manual(values =colsdm) +
  ylim(1,2) +
  xlab("")+
  ylab("99th percentile / Median")

plot(plot2)
filenamesdm <- paste("../output/png/return/",yearname,regionname,varname,"return_SDM.png")
ggsave(plot2,filename=filenamesdm, dpi=250,width=7.5, height=6)

}
}


  
#---end return period


# Figure2b
#----- #4 PPOP_UNSH & PCAL: Bars across SSPs & regions ----
var.in <- c("PCAL")
var.in <- c("PPOP_UNSH")
var.in <- c("PCAL","PPOP_UNSH")

for(k in 1:length(var.in)){
  varname <- var.in[k]
  
  
if (varname=="PCAL"){
  labelname <- c("calorie intake [kcal/person/day]")
} else if (varname=="PPOP_UNSH"){
  labelname <- c("population at risk of hunger [billion]")
} else if (varname=="YEXO"){
  labelname <- c("yield [tonne/ha]")
} else {
  labelname <- c("aaaaaaaa")
}
#colnames <- c("VAR","Year","Region","SSP","GCM","RCP","CO2fert","SYR","N","Crop","Value")

    y <- pdata[pdata$VAR==varname & pdata$Year==yearname  & pdata$RCP=="RCP85" & pdata$CO2fert=="noco2" & pdata$N=="N1" & pdata$Crop=="dum" & pdata$Region!="WLD", c(-1,-2,-6,-7,-10)]

	names(y) <- c("Region","SSP","GCM","SYR","N","Value")	#Rename the column

	minv <- min(y[6])
	maxv <- max(y[6])

	interact<-interaction(y$SSP, y$Region, sep=" : ")
	
#	bw <- diff(range(y$Value))/50
#	mu <- ddply(y, "RCP", summarise, grp.mean=mean(Value))
#	mu1 <- ddply(ypre, "RCP", summarise, grp1.mean=mean(Value))


  y0 <- ddply(y,c("interact"),transform,v0=min(Value))
  y1 <- ddply(y0,c("interact"),transform,v50=median(Value))
  y2 <- ddply(y1,c("interact"),transform,v100=max(Value))
	z <- regionlabel(y2)
#	z <- regionlabelagg(y2)

gbox <- ggplot() + 
      geom_bar(data=z, aes(x=as.factor(Region),y=v50), position=position_dodge(), stat="identity", colour='black', fill='grey60') +
      geom_errorbar(data=z, aes(x=as.factor(Region), ymin=v0, ymax=v100), width=.2,position=position_dodge(.9)) +
#    geom_bar(data=z, aes(x=as.factor(Region),y=v50, fill=SSP), position=position_dodge(), stat="identity", colour='black') +
#    geom_errorbar(data=z, aes(x=as.factor(Region), ymin=v0, ymax=v100, fill=SSP), width=.2,position=position_dodge(.9)) +
#    geom_boxplot(data=z, aes(x=Region,y=Value, color=interact, ymin = v0, lower = v25, middle = v50, upper = v75, ymax = v100), fill="grey80", stat="identity") +
#  	ylim(minv-abs(minv)*0.01, maxv+abs(maxv)*0.01) +
#	ggtitle(mainlabel)
	ylab(labelname) + xlab("") +	MyThemeLine
#  geom_point(data=z, aes(x=Region, y=NoCC),size=10,shape=73,colour="black")

  gbox1 <- gbox + coord_flip()
  plot(gbox1)

  gbox2 <- gbox1 + facet_wrap(~SSP,ncol=2) 
  #+ scale_color_manual(values =col)
  plot(gbox2)
  
  filenamebox <- paste('../output/png/', varname ,'/',yearname,'barsspagg.png', sep = '')
  ggsave(gbox2,filename=filenamebox, dpi=250,width=10, height=4)

}

#----End --- PPOP_UNSH & PCAL: Bars across SSPs & regions ----

#-----  #5 YIELD(YEXO) & price(XPRP): Boxplot across regions ----


var.in <- c("XPRP")
var.in <- c("YEXO")

for(k in 1:length(var.in)){
  varname <- var.in[k]
  
  if (varname=="PCAL"){
    labelname <- c("calorie intake [kcal/person/day]")
  } else if (varname=="PPOP_UNSH"){
    labelname <- c("population at risk of hunger [billion]")
  } else if (varname=="YEXO"){
    labelname <- c("yield [tonne/ha]")
  } else if (varname=="XPRP"){
    labelname <- c("production price [2005=1]")
  } else {
    labelname <- c("aaaaaaaa")
  }
  
#  filenamebox2 <- paste('../output/png/', varname ,'/',yearname,'_limited.png', sep = '')
  filenamebox2 <- paste('../output/png/', varname ,'/',yearname,'.png', sep = '')
  
  #colnames <- c("VAR","Year","Region","SSP","GCM","RCP","CO2fert","SYR","N","Crop","Value","NoCC")
  
  # Limit the num of regions
  #pyld <- pdata[pdata$VAR=="YEXO" & pdata$Year==yearname  & pdata$SSP=="SSP2" & pdata$Region!="WLD" & 
  #               pdata$Region!="XOC" & pdata$Region!="XNF" & pdata$Region!="XER" & pdata$Region!="XE25" & 
  #               pdata$Region!="USA" & pdata$Region!="TUR" & pdata$Region!="JPN" & pdata$Region!="CIS" & 
  #               pdata$Region!="CAN"& pdata$Region!="BRA" & pdata$Region!="CHN", c(-1,-2,-4)]
  
  # Limit the num of regions
  pyld <- pdata[pdata$VAR==varname & pdata$Year==yearname  & pdata$SSP=="SSP2" & pdata$Region!="WLD", c(-1,-2,-4)]
  
  
  names(pyld) <- c("Region","GCM","RCP","CO2fert","SYR","N","Crop","Value","NoCC")	#Rename the column
  
  minv <- min(pyld[8],pyld[9])
  maxv <- max(pyld[8],pyld[9])
#  maxv <- c(20)
  
  #	bw <- diff(range(y$Value))/50
  #	mu <- ddply(y, "RCP", summarise, grp.mean=mean(Value))
  #	mu1 <- ddply(ypre, "RCP", summarise, grp1.mean=mean(Value))
  
  interact<-interaction(pyld$RCP, pyld$CO2fert, sep=" : ")
  
  pyld0 <- ddply(pyld,c("Region","interact","Crop"),transform,v0=min(Value))
  pyld1 <- ddply(pyld0,c("Region","interact","Crop"),transform,v25=quantile(Value,0.25))
  pyld2 <- ddply(pyld1,c("Region","interact","Crop"),transform,v50=median(Value))
  pyld3 <- ddply(pyld2,c("Region","interact","Crop"),transform,v75=quantile(Value,0.75))
  pyld4 <- ddply(pyld3,c("Region","interact","Crop"),transform,v100=max(Value))
  pyld5 <- regionlabel(pyld4)
  pyldz <- croplabel(pyld5)
  
  
  gyldbox <- ggplot() + 
    geom_boxplot(data=pyldz, aes(x=Crop, y=Value, color=interact, ymin = v0, lower = v25, middle = v50, upper = v75, ymax = v100), fill="grey80", stat="identity") +
    ylim(minv-abs(minv)*0.01, maxv+abs(maxv)*0.01) + ggtitle("") +
    ylab(labelname) + xlab("") +	MyThemeLine +
    geom_point(data=pyldz, aes(x=Crop, y=NoCC),size=6,shape=73,colour="black")
  
  gyldbox1 <- gyldbox + coord_flip()
  gyldbox2 <- gyldbox1 + facet_wrap(~Region, ncol=4)
  
  plot(gyldbox2)
  
  ggsave(gyldbox2,filename=filenamebox2, dpi=250,width=7.5, height=7.5)
  
}



#----- END --- YIELD : Boxplot across regions ----






